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Abstract 

We present an algorithm for computing some previously introduced concepts of rel- 
(~| ative velocity of a distant test particle with respect to an observer: kinematic, Fermi, 

spectroscopic and astrometric relative velocities. The algorithm is focused on finding the 
'"5 Fermi and observational coordinates of the test particle with respect to the observer, and 

it can be extended to non-convex normal neighborhoods. We make some convergence 
T-H tests and study some fundamental examples in Schwarzschild and Kerr spacetimes. Fi- 

nally, we show an alternative method for computing the Fermi and astrometric relative 
^ velocities. 

cr 

^ 1 Introduction 

In Newtonian mechanics, the concept of "relative velocity" of a test particle with respect 
" ' to an observer is unambiguous and fundamental. Nevertheless, in general relativity it is 

only well defined when the observer and the test particle are in the same event. In order to 
generalize this concept for distant test particles, a definition of relative velocity based on a 
particular coordinate system can be introduced, and it could be interesting in some cases, like 
the comoving coordinates in FLRW spacetimes for static observers (see [^). Moreover, some 
notions of relative velocity of a distant test particle were introduced by the lAU using adapted 
reference systems in the case of objects in the neighborhood of the solar system (see [2][3]). 
However, a universal concept of relative velocity independent from any coordinate system is 
' . I basic and thereby, some authors have proposed geometric definitions without any coordinate- 

J> dependence (see [4[|6]). In this way, four different intrinsic geometric definitions of relative 

k>( velocity of a distant test particle with respect to a single observer were introduced in [?]. 

These definitions are strongly associated with the concept of simultaneity: kinematic and 
Fermi in the framework of "spacelike simultaneity" , spectroscopic and astrometric in the 
framework of "lightlike simultaneity" . 

These four concepts of relative velocity each have full physical sense, and have proved to be 
useful in the study and interpretation of properties of particular spacetimes (see [7}|l2]). But, 
in most cases, the computations are analytically very complex and this makes the theoretical 
study much harder. So, we present in this paper a numerical algorithm that enables the 
computation of these relative velocities. 

This paper is organized as follows. In Section[2]we present the framework, establishing the 



notation and defining some necessary concepts, introducing in Section 2.1 the four geometric 
concepts of relative velocity and extending the original definitions to non-convex normal 
neighborhoods. In Section |3] we develop the algorithm focused on finding the Fermi and 
observational coordinates of the test particle with respect to the observer by means of finding 
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a geodesic (with certain characteristics) from the observer to the test particle, and we also 
make a discussion about the convergence of the method in Remark |3.4[ In Section |4] we 
give some fundamental examples in Schwarzschild and Kerr spacetimes, showing the rate of 
convergence and the effectiveness of the algorithm. Finally, in Appendix |A] we present an 
alternative method for computing the Fermi and astrometric relative velocities. 

2 Definitions and notation 

We work in a Lorentzian spacetime manifold {Ai,g), with c = 1 and V the Levi-Civita 
connection, using the "mostly plus" signature convention (—,+,+,+). Given two events 
p, q, and a segment curve "0 that joins p and q, the parallel transport from p to q along 
■0 is denoted by r^. Given a curve /3 : I ^ Ai with / CM, the image /?/ (a subset in 
A4) is identified with /3. Vector fields are denoted by uppercase letters and vectors (defined 
at a single point) are denoted by lowercase letters. If u is a vector, then denotes the 
orthogonal space of u. The projection of a vector v onto is the projection parallel to u, 
i.e. V — u- Moreover, if a; is a spacelike vector, then ||x|| :— g (a;,x)^^^ is the modulus 
of X. Given a vector field X, the unique vector of X in TpM is denoted by Xp. 

In general, we say that a timelike world line /3 is an observer (or a test particle); never- 
theless, we say that a future-pointing timelike unit vector u in TpA4 is also an observer at p, 
identifying the observer with its 4- velocity. 

A light ray is a lightlike (null) geodesic A. A light ray from q to p is a, light ray A such 
that q^p € X and p is in the causal future of q. 

We are going to consider two kinds of intrinsic simultaneity: spacelike and lightlike. Given 
an observer u at p, the events simultaneous with u form the corresponding simultaneity 
submanifold: 

• Spacelike simultaneity: the Fermi surface Lp .^ (also known as Landau submanifold) 
is given by all the geodesies starting from p and orthogonal to u. In terms of the 
exponential mapj^] on TpM. , it is given by expp u-^ . 

• Lightlike simultaneity: the past-pointing horismos submanifold Ep is given by all the 
light rays arriving at p, i.e. it is given by exp^ Cp where C~ is the past-pointing light 
cone in TpM composed by all the past-pointing lightlike vectors of TpM. 

2.1 Geometrically defined relative velocities 

Four different definitions of relative velocity of a test particle with respect to an observer 
were introduced in [t], working in a convex normal neighborhood, where given two different 
events there exists a unique geodesic joining them. Now, we are going to work in a general 
spacetime M, not necessarily a convex normal neighborhood, and we are going to extend the 
definitions of [7j to this new setting, where two events could be joined by more than one (or 
none) geodesic and the simultaneity submanifolds could present self-intersections. 

Throughout the paper, we consider an observer /3 and a test particle /3' (parameterized 
by their proper times) with 4-velocities U and U' respectively. Moreover, we consider an 
event p oi (3 with 4- velocity u := Up. 

Given a vector s ^ such that expp s e /?', the corresponding kinematic relative velocity 
of 13' with respect to u is defined by the vector 

Vkin -r^ \'^tp'^'s - (1) 

-9 [T-Zpu'^,uj 

* Given V £ TpAi, cxp^ v := 7^(1) where 7^ is the geodesic starting at p with initial tangent vector v. 
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Spacelike simultaneity Lightlike simultaneity 




Figure 1: Scheme of the elements involved in the study of relative velocities of with respect 
to u in general (not necessarily in a convex normal neighborhood). Left: given s S such 
that expp s € /?', we define ^p (a) := exp^ as for < a < 1, qg := "0(1)7 and u'^ is the 4- velocity 
of at Qs- Right: given w € C~ such that exp^ w € /?', we define Sobs ■— w + g{u,w)u, 
A (a) := exp aw for < a < 1, := A(l), and u'^ is the 4- velocity of /?' at g£. 



where qs '■— expp s is the event of /3' at which the relative velocity is measured, u'^ :— U^^ is 
the 4- velocity of /?' at q^, and r^j, is the parallel transport from qs to p along the geodesic 
segment given by ip {a) :— exp^ as for < a < 1 (see Figure [IJ left). In this case, s is 
a relative position of with respect to u, and there is a different Wkin for each different s 
satisfying s € u^, exp^ s £ /3'. Note that if we work in a convex normal neighborhood then 
s is unique. 

Analogously, we can define another concept of relative velocity also introduced in [4| : given 
a vector w e C~ such that expp w € /?', the corresponding spectroscopic relative velocity of 
P' with respect to u is defined by the vector 



-9 (Tqepu'e^u) 



X 



.u'e - u. (2) 



where qi :— exppW is the event of /?' at which the relative velocity is measured, :— U^^ 
is the 4- velocity of /?' at qi, and r^^p is the parallel transport from q^ to p along the light 
ray segment given by A (a) := exp^ aw for < a < 1 (see Figure [Tj right). In this case, the 
projection of w onto u-^ given by Sobs '■— w + g{u,w)u is the corresponding observed relative 
position of j5' with respect to u, and there is a one-to-one correspondence between w and Sobs 
(note that w = Sobs — ||sobs||u). So, there is a different Vspec for each different w satisfying 
w e C~, exppW g /3', and this means that if there is gravitational lensing then each image 
of the observed object has a different spectroscopic relative velocity. 

Remark 2.1 The spectroscopic relative velocity is specially useful because the frequency 
shift can be deduced from it (see (?]): let A be a light ray from qi to p and let u, be two 
observers at p, qi respectively; then 



i 1,^1 Sobs 



(3) 



where v' are the frequencies of A observed by u, respectively, Wgpoc is the corresponding 
spectroscopic relative velocity given by ([2]), and Sobs is the corresponding observed relative 
position of v!n with respect to u. 
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In general, the existence of s or w (and consequently Sobs) is not assured. Nevertheless, 
if they exist for each event of the observer /3, we can construct (difFerentiable) vector fields S 
and S'obs defined on j3, representing a relative position and an observed relative position of f3' 
with respect to /?, respectively. But, of course, not all choices of s and w lead to differentiable 
S and 5obsj only those that change differentiably when varying the event of the observer. 
From S and S'obs, we can construct the vector fields Vkin and Vspcc defined on /?, representing 
a kinematic and a spectroscopic relative velocity of /?' with respect to j3, respectively. 

Given a vector field S defined on f3 and representing a relative position of /?' with respect 
to /3 (i.e. such that Sp € and exp^ Sp € /?' for all p G /3), the corresponding Fermi relative 
velocity of /?' with respect to /3 is the vector field 

V^Formi + g {VuS, U)U^VuS-g (5, VuU) U, (4) 

defined on /3. 

Analogously, given a vector field S'obs defined on /3 and representing an observed relative 
position of /?' with respect to /3 (i.e. such that Sobs p G andexpp(Sobs p~||Sobs pIlC^p) G P' 
for all p Cz (3), the corresponding astrometric relative velocity of /?' with respect to j3 is the 
vector field 

Kst := V[/ S'obs + g (V[/ S'obs, U)U ^ Vt/Sobs - 9 (Sobs, V[/C/) C/, (5) 

defined on /3. 

If we work in a convex normal neighborhood, then it is assured that there exists a unique 
S and Sobs, and hence there exists a unique Vkin, Kspcc, V^Fcrmi, and Vast (see |7j). But this is 
not true in general and so different vector fields S and Sobs define different vector fields T4in, 

V^poc, ^^Fcrmi, and Vast- 

In order to complete the notation that we are going to use, we define the vectors WFcimi '■= 
Vpermip and Uast := V^stp] morcovcr, throughout the paper we are going to denote s := Sp, 
Sobs S'obs pi 'I'kin := Vcinp, and ■^spec := Vspecp as we have already done in this section. 



3 Algorithm for computing geometric relative velocities 

First, we are going to suppose that we work in a convex normal neighborhood, and later, in 



Section 3.4 we will extend the discussion to non-convex normal neighborhoods. So, working 
in a convex normal neighborhood implies that given an event p € P with 4-velocity u, there 
exists a unique event Qs G /?' such that the unique geodesic ip that joins p and is in Lp „ 
(i.e. it is orthogonal to u at p); on the other hand, there exists a unique event qi £ /3' such 
that the unique geodesic A that joins p and qi is in Ep (i.e. it is a light ray arriving at 
p). In this case, the main difficulty is to find the geodesies V' (spacelike simultaneity) and A 
(lightlike simultaneity), taking into account that the events qs and qi are also unknown (see 
Figure [I]). This is equivalent to finding the relative positions s and Sobs, that are determined 
by the Fermi and observational (or optical) coordinates (see 13-17]). Once we have found 



them, the corresponding relative velocities are easy to compute by means of their definitions 



(see Section 2.1 1 



If we can not find theoretically the Fermi and observational coordinates, then we can not 
find the geodesies ip and A of Figure [T] Of course, there is a brute-force method for finding 
these geodesies that consists on launching a lot of geodesies from the observer at p in different 
directions (orthogonal to u in the spacelike case, and past-pointing lightlike directions in the 
fightfike case), trying to reach the test particle /3' . But obviously, this method spends a lot 
of computation time, and solving the geodesic equations with an acceptable accuracy is not 
as fast as we desire in some metrics. 

So, we are going to propose a more efficient method based on an iterative correction 
algorithm with a Newton- Raphson structure. But first, we need to make more precise the 
concept of "nearness" between curves. 
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3.1 The concept of "nearness" 

There is no global concept of distance in pseudo-Riemannian manifolds because the metric g 
is degenerate. Nevertheless, given a 3-dimensional spacelike foliation in M. we have that the 
induced metric g is Riemannian. Hence, we are going to suppose that we use a coordinate 
system {a;°, a;^, a:^, x'^} such that x^, x^, are spacelike coordinates and a;° (also denoted as 
t) is a timelike coordinate, referred as coordinate time] then, t = constant defines the leaves 
of the desired spacelike foliation, and the covariant coefficients of the induced Riemannian 
metric g are given by g.^- = gij , where g^^ are the covariant coefficients of the general metric 
g in the above coordinate system (Latin indices run over 1,2,3, and Greek indices run over 
0, 1, 2, 3). Namely, given two vectors f , w in the same tangent space of a leaf, we have that 
g{v,w) = gijV^w^ , and ||u|| — g{v,vy^^ . 

Therefore, given any event p in a leaf, there exists a local concept of spatial distance for 
events q in the same leaf: d(p,q) :— ||w||, where v is the vector which "joins" p and q, i.e. 
such that exppW = q where exp^ is the induced exponential map at p. If d{p,q) « 0, the 
tangent spaces at p and q can be identified by means of the coordinate system, and we can 
consider an affine structural around p, having exppV ~ p + v, i.e. « — p*. So, we will 
say that two events p,q with the same coordinate time are close if \\q — p\\ is considered to 
be small, where q — p is the vector in the tangent space of p with coordinates — p^- 

This concept of nearness can be also applied to curves: we will say that two curves c, c' 
are close if there exist two events p d c and g G c' with the same coordinate time such that 
p and q are close (i.e. \\q — p\\ is small). 

3.2 Spacelike simultaneity 

Let {t ^} be a coordinate system such that x^,x'^,x^ are spacelike coordinates 

and t is a timelike coordinate. In the framework of spacelike simultaneity and taking into 
account the previous concept of nearness, we choose an initial vector sq G w^, such that 
the geodesic -00 starting from p with initial tangent vector sq is sufRciently close to the test 
particle /3'; i.e. there exist events ggeo (in the geodesic) and (?part (in the test particle), both 
with the same coordinate time, such that ||(7part — 9gco|| is considered to be small. 

Remark 3.1 In practice, we need a sub-algorithm for estimating ^goo and qp^rt- given the 
initial geodesic "001 we can compute the spatial distance between events of ipo ^^nd the corre- 
sponding events of /?' with the same coordinate time; supposing that qgco (in V'o) and (/part (in 
/?') are the events that minimize this distance, we can use a bisection method for estimating 
them with low computational cost. 

Since we are going to work near the event ggeo, we can identify all the tangent spaces by 
means of the coordinate system and provide an affine structure in the vicinity of Qgeo) where 
exp^^^^ V ~ qgco + v. Hence, the Fermi surface Lp^u can be linearly approximated by the affine 
hyperplane ^goo + Tq^^^Lp_u and the intersection event 5s between and q^co + Tq^^^Lp^u is 
approximated by 

qs ■= ggco + h, (6) 

where h is the projection of gpart — 9gco onto Tq^^^Lp^u parallel to Upart: with Upart the 4- 
velocity of (3' at gpart (see Figure [2] with n — 0). 

For finding h we need first to find Tq^^^Lp^u, and for this purpose we can use the Jacobian 
matrix Jj^ := d^, exp^ sq where 9^ is the partial derivative with respect to the ^-th coordinate 
(note that the derivatives of exp^ are easy to estimate numerically). So, {ci, e2, 63} is a basis 
of Tq^^^Lp^u, where ef := J^e^ and {ei, 62, 63} is a basis of u-^, e.g. the projections onto u-^ 

t Given an affine structure, we can subtract points to get vectors, or add a vector to a point to get another 
point. 
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Figure 2: Diagrams for Qs (left) and the estimations of qs (right) when ggeo and gpart (that 
are events with the same coordinate time) are close. In this case, all the tangent spaces are 
identified by means of the coordinate system and provide an afRne structure around qgeo- 
Right: the Fermi surface Lp_„ is approximated by the affine hyperplane ggeo + Tg^^^Lp^u- 



of the spacelike vectors 



dx^ I p 



(0,1,0,0), 



(0,0,1,0) and 



d I — 



= (0,0,0,1). Then, 



supposing that the Fermi surface Lp^u is spacelike at qgeo, we have that {ei, 62, 63, Mpjj,.t} is 
a basis of Tg^^^A^, and hence h — a^Ci, where a* are the spatial coordinates of gpart — qgeo in 
the above basis. 



Remark 3.2 Given q £ Lp ^, in 18 Proposition 3] it is proved that, in some cases, TqLp ^ = 
(jpqU)^ and consequently, Lp_„ is spacelike at for example, in the case of stationary ob- 
servers in the Schwarzschild spacetime (see Examples 4.1 and 4.3). Nevertheless, TqLp u is 



not (Tpqu) in general, as it occurs in the analogous case of stationary observers in the Kerr 
spacetime (see Examples 
remains spacelike. 



4.2 



and 4.4), although it can be checked that the Fermi surface L 



p,u 



The next step is finding qg by means of a Newton-Raphson method: we find qg by 
and we redefine gpart as the event of f3' with the same coordinate time as 5s; then, by Ml) 
again, we find another qg, repeating this process until qs approximates 5s with the desired 
accuracy. The convergence of this method is assured by the next proposition. 



Proposition 3.1 If the Fermi surface Lp ^ is spacelike at qgeo and the acceleration of the test 
particle (3' is bounded, then the Newton-Raphson method described above produces a sequence 
of events {^snlngN that converges to qs with quadratic order for a sufficiently close events 

gpart and qgeo- 



Proof. We are going to denote 5s n as g„ for the sake of simplicity. Working in coordinates 
{t, , , x^} , let (3' be parameterized by the coordinate time t and let to '-— ^part: '-— 9n-i 
for n > 1. 

Given n e N, the event q„ is the intersection of the afhne hyperplane qgeo + Tq^^^Lp .^ 
and the affine line /3'{tn) + (/3'(i„)), where the overdot denotes derivation with respect to t 
and "( )" denotes the span. If t + OiX^ = 6 is the cartesian equation of qgeo + Tq^^^Lp .^ (we 
can suppose that the coefficient of t is 1 because Lp^u is not timelike), then, by algebraic 
manipulations, we have that the coordinate time of q„ is given by 



(tjiXj-^ 



x^ 



1 



(7) 



where xjj :— f3'^{tn) and i*^ := /3'*(i„) for i = 1,2,3. It can be proved that 1 + a^i^ ^ 
because /?' is timelike and Lp^u is spacelike at qgeo, and so the intersection event q„ always 
exists. 
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On the other hand, if qs — /3'{i), from the Taylor expansion of order 1 of (3'^{t) at i„ we 
have that 

qi ^ xl, + xl,e„ + R{, (8) 

where e„ i— t„ and R\ := 5/3'n?n)£n with between i and t„ for i — 1,2,3. Moreover, 
since q^ G qgeo + Tq^^^Lp ^, we have that 

i + a.qi = b. (9) 

Substituting Q in (|9]) and taking into account ([t]), we get 



,2 



1 + a.,xl 2(1 + a^ij,) 



Assuming that the acceleration of (3' is bounded and |e„| is sufficiently small (this can be 
achieved if qp^rt and ggco are sufficiently close), the result holds. 

■ 

Once we obtain qs, we have to estimate the vector s^such that the geodesic starting from 
p with initial tangent vector 's arrives at qs, i.e. s := exp"-'^ 5s. Note that ^ might not be 
exactly orthogonal to u because qgco + Tq^^^Lp^u does not coincide in general with (i.e. 
qs is an estimation of qs and they do not coincide in general). Re-scaling sq in order to 
verify expp sq = ggoo (for convenience) and working in coordinates, we can make the linear 
estimation 

qs^^q^co + Jl^-is'-s^o), (10) 

where = di, exp^ sq was previously computed. Hence, solving the linear system given by 



( 10 1 we have 



s'(:^s^, + iJ-X-&'~q',eo)^s'', (11) 

where (J~^)J^ are the coefficients of the inverse of the Jacobian matrix J^. 

Finally, since Si might not be exactly orthogonal to u (as it happens with s), we have to 
redefine si projecting it onto u-^, obtaining in this way an estimation of the desired vector s. 

Remark 3.3 A higher order approximation based on the Taylor expansion of exp can be 



applied instead of the linear estimation (10 1. For example, the quadratic approximation 



ql" ~ <eo + exp^ so ■ (s - s^) + ^9,, exp^ sq • {s" ~ s^) ' " - s^)- (12) 



In this case, the coordinates of s can be solved from the system of quadratic equations (12), 



obtaining an expression equivalent to (111 but involving a second order approximation. If 
there is more than one solution for s, we have to take the solution such that exp^ 's is "closest" 
to qs, e.g. that minimizes X)/l=i I exp^ s"— q^\. However, in practice, usually it would suffice 
to take the closest (in a coordinate sense) solution to sq. 

Let tpi be the geodesic with initial direction si. Probably, this new geodesic V'l does not 
intersect (3' because we have made estimations, but it should be closer to /3' than the initial 



geodesic ipo (see Remark 3.4). So, applying this method again to the new geodesic, we get 
other values for (?gco, ^part, and qs, from which we obtain S2 and a geodesic ip2 (with initial 
direction S2) that should be closer to f3' than ipi. We must repeat this process, obtaining 
geodesies ipn (with initial directions s„), until we reach the desired accuracy (see Figure [2]) 
or we exceed a certain number of iterations. 

Summing up, given an observer's event p = /3{t) with 4-velocity u, the algorithm for 
finding the relative position s of the test particle /3' is as follows: 

1. Choose an initial vector sg £ such that the geodesic i/'o starting from p with initial 



direction sq is sufficiently close to the test particle /?' (in the sense of Section 3.1 ). Set 
n = 0. 
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Find ggoo G V'n and (/part G P' with the same coordinate time that minimizes the 



distance ||(7part — 9gco|| (see Remark 3.1 1. For n > 1 we must check that this distance is 
lesser than the corresponding distance computed for n — 1; if it does not hold, then we 
should stop the algorithm because it could be a divergence symptom. 

Find Qs by means of a Newton- Raphson method: we find by (|6|, and we redefine 
9part as the event of /?' with the same coordinate time as qs] then, by ^ again, we find 
another gg, repeating this process until qs approximates qs with the desired accuracy. 
A quadratic convergence is assured by Proposition |3.1[ 



4. Define s'^^^^ := sj^ + ( J ■ (^q^ — 9^0) (it is a linear approximation of s :— exp ^ q^ 



see (11)), where s„ is re-scaled in order to hold exp^ s„ = ggco. Alternatively, for 



a second order approximation, we can define Sn+i by solving a system of quadratic 



equations (see Remark 3.3 ) 



5. Redefine Sn+i as its projection onto u^: s„+i — Sn+i + g{u, s„+i)w. 

6. Set n = n + 1. Repeat the process (steps 2, 3, 4, 5) with the geodesic ipn starting from 
p with initial direction s„, until we reach the desired accuracy. Otherwise, we should 
stop the algorithm if we arrive at a predetermined maximum number of iterations. 

If the desired accuracy has been achieved, we can apply this algorithm again for another 
event P{t + At) of the observer. Then, the new initial vector sq should be chosen as the vector 
with the same coordinates as the corresponding final vector s„ computed for the previous 
event f3{t). In this case, choosing a sufficiently small time step At assures convergence (see 



Remark 3.4), and "differentiability" of S (although we obtain a discrete vector field) in the 



case of non-convex normal neighborhoods (see Section 3.4 1 



Remark 3.4 With respect to the convergence, in practice, this algorithm has been tested in 
many spacetimes with different test particles and observers, and it numerically converges in 
all of them with quadratic order. This accords with the fact that, if we compute qs exactly 
(in the 3rd step) and define s„+i :— s (in the 4th step), the algorithm has a Newton- Raphson 
structure. 

But, in general, it is not assured that ipn+i is closer to /3' than ipn- If this property does 
not hold, then convergence is not guaranteed. Determining theoretical conditions for assuring 
convergence in a general case is a hard open problem. For this purpose. Proposition 3.1 is 
the first step. 

Nevertheless, if the algorithm converges at the observer's event p = /3(t), then, applying 
differentiability arguments, it is assured that there exists a sufficiently small time step Ai > 
such that the algorithm also converges at f3{t + At) taking the new initial vector sq as the 
vector with the same coordinates as the corresponding final vector s„ computed for the 
previous event p. 

3.3 Lightlike simultaneity 

Analogously, in the framework of lightlike simultaneity, the initial vector wq £ TpAi must 
be liglitlike and past-pointing, and the initial geodesic is named Aq. Then, we estimate the 
intersection of the test particle with the affine hyperplane qg^o '^Tq^^^E~ using an expression 
analogous to ([6]): 

qt ■= ggco + (13) 

where h is the projection of (/part ~ 'Zgoo onto Tq^^^E~ parallel to Wpart? with Wpart t he 4- velocity 
of /?' at (/part (see Figure |3] with n = 0). Since is hghthke at q (see Remark 3.5 below), 
we have that h is always well-defined. 
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Figure 3: Diagrams for qi (left) and the estimations of qi (right) when qg^o and (/part (that 
are events with the same coordinate time) are close. In this case, all the tangent spaces are 
identified by means of the coordinate system and provide an affine structure around qgeo- 
Right: the past-pointing horismos submanifold E~ is approximated by the affine hyperplane 
<?gco + Tq^^^E~ . Note that Tq^^^E~ = w^^, where Wgco is the tangent vector of A„ at ggoo- 



Remark 3.5 Given g G i?p , in 18 Proposition 3] it is proved that TqE^^ = (rp^exp^ ^ q)^ 
and hence, the past-pointing horismos submanifold is always lightlike; in fact, we can write 



TqEp 



w where w is the tangent vector of the light ray from q to p. So, 



?(^gcoj (/part 9gco) ^/ 



9{wg 



t) 



part ' 



(14) 



where Wg^o is the tangent vector of Ag at ^gco- Note that g{w^ 
lightlike and Up^rt timelike. 



gco I ■'^part 



) 7^ because Wg^o is 



Next, we have to find q^ applying a Newton- Raphson method analogously to the spacelike 



case. It can be proved analogously to Proposition 3.1 that this part of the algorithm has a 
quadratic order of convergence, but we only need the hypothesis on the bounded acceleration 
because El^ is always lightlike (see Remark 



3.51 



After this, we have to estimate the vector w := exp^ ^ qg. For this purpose, we re-scale Wq 
in order to verify exppWo = ggeo and, working in coordinates analogously to (111, we obtain 



Wn 



J- 



w 



(15) 



where, in this case, (j are the coefiicients of the inverse matrix of (9j, exp^(wo). 

A second order approximation can be deduced analogously to Remark |3.3| replacing, in ex- 
pression (12), Qs, s", and Sq with qi, w, and respectively. 



Finally, since wi might not be exactly lightlike (as it happens with w) , we have to redefine 
Wi projecting it onto C~ , obtaining in this way an estimation of the desired vector w. 
For example, we can redefine it as the past-pointing lightlike vector with the same spatial 
coordinates as the original wi (i.e. a projection parallel to ■§i\p)- 

The steps of the algorithm in the framework of lightlike simultaneity are analogous to 
those exposed at the end of Section [3^2] in the framework of spacelike simultaneity, replacing 



V'ra, ?s, and qs with w„, , A„, and qi respectively. In the 3rd step, we can 



compute q^ by means of expression (14). Moreover, in the 5th step, we have to project Wn+i 



onto Cp as it is explained in the above paragraph. 
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3.4 Non-convex normal neighborhoods 



Working in a convex normal neighborhood, there is no problem in the determination of the 
events qgco (in the geodesic ipn or A„) and (/part (in the test particle /?') introduced in the 
beginning of Section 1X2] because they globally minimizes the distance between events of the 



geodesic and the test particle in surfaces of constant coordinate time (see Remark 3.1 ). But 
if we work in a non-convex normal neighborhood, then there could be different (or none) 
possibilities of relative position of the test particle with respect to the same event p of the 
observer and, in this case, each possibility of relative position drives to a different local 
minimum of this distance. Hence, for determining a "suitable" pair Qgco, 9part we have to 
search the local minimum that corresponds to a "suitable" relative position, according to the 
previously computed relative positions. 

Let us explain it in more detail: suppose that we have previously applied the algorithm 
in pq := (3{t) and we have obtained a relative position Sp„. Then, applying the algorithm in 
Pi :— P{t + Ai), we say that a relative position s^^ is suitable (according to Sp^) if Sp^ and 
Sp^ are vectors of a (differentiable) vector field S of relative positions. If the time step Ai is 
sufficiently small, then it is assured that there exists at least one suitable Sp^ . If there are 
several relative positions, a suitable one must be close (in a coordinate sense) to Sp„ , and so 
it can be chosen as a relative position that, for example, minimizes the coordinate distance 
S^=o l^pi ~ *pol' there are several suitable relative positions, then any of them are valid. 
In this case, the output can be controlled adding some desired restrictions. For example, in 
the case of an equatorial circular geodesic test particle and an equatorial stationary observer 
in Schwarzschild spacetime (where there is gravitational lensing, see Example 4.3), it could 
be desirable that all the involved geodesies ('0n and A„) were also equatorial; we can impose 
this condition and so the relative positions S and Sohs have zero ^-component. 

In practice, a suitable pair ggco, ^Zpart for pi satisfies that the sum of the coordinate dis- 
tances between them and the corresponding pair for po is small, given a sufficiently small 
time step At. 

Taking all this into account, the algorithm presented in Section [3] returns one possibility 
of a discretized version of a relative position vector field S or S'obs, provided that there exists 
a relative position in the first event of the observer in which we apply the algorithm and we 
use a sufficiently small time step in the observer. From this output we obtain one discretized 
version of the (differentiable) vector fields 14;in, Vpermi, Kspec, or Vast- 



4 Examples 

The algorithm has been tested in several spacetimes with different observers and test particles, 
obtaining very good results in computation time and accuracy. We present here the most 
representative examples in Schwarzschild and Kerr spacetimes, using a test particle with 
equatorial geodesic orbit and stationary observers. 

To sum up, given an observer and a test particle, our objective is to find s (spacelike 
simultaneity) or w (lightlike simultaneity) at a given event p of the observer (see Figure [T]). 
To do this, we need an initial vector so or wq that is supposed to be close (in a coordinate 
sense) to s or w respectively; but first, we are going to show by means of Examples |4.1| and 
|4.2| the rate of convergence of this method using an initial vector not necessarily close to the 
objective vector. Moreover, we are going to apply the second order approximation proposed 
in Remark |3.3| and compare with the usual linear method. 

Example 4.1 The Schwarzschild metric in spherical coordinates {t,r,9,if} is given by the 
line element 

ds'- ^ - ("l - ^"j di^ + ("l _ ^I!^") ' + r2 (d02 + gij^2 Q^^2^ ^ (ig) 
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Spacelike simultaneity 



Lightlike simultaneity 




Figure 4: Diagrams in the a;y-plane of the elements involved in Example 4.1 The sta- 



tionary observer launches an initial geodesic from p — (0, 8, 7r/2, 0) with initial direction 
So = (0,0,0,1) (left) and wq = (—^,0,0,1) (right). It is shown how the successive iter- 
ations of the algorithm return geodesies that are getting closer to the desired intersection 
point with the test particle, or (see Figure [T]). 



where the parameter m is interpreted as the mass of the gravitating object, r > 2m is the 
radial coordinate, and Q < 9 < tt. From now on we are going to suppose that the coordinates 
hold these restrictions and m = 1. In the framework of this coordinate system, a stationary 
observer is an observer with constant spatial coordinates. Note that stationary observers are 
not geodesic, but they are useful in the description and interpretation of the Schwarzschild 
spacetime. 

Given a stationary observer at j'q = 8, = 7''/2, <y3o — 0, and a test particle with 
equatorial circular geodesic orbit with ri — A, 9i = 7r/2, ipi — n jl at t = 0, we are going 
to check the algorithm for computing s and w at p = (0, 8, 7r/2, 0) using an initial spacelike 
direction sq = (0,0,0, 1) (orthogonal to the 4-velocity of the observer at p), and using an 
initial past-pointing hghtlike direction = (— ^tq, 0, 0, 1), respectively. Note that these 
initial directions are not too close of the objective vectors s and but despite this, the 



method works well. Moreover, applying the second order approximation of Remark 3.3 gives 
better results, but it doubles the computation time because we have to solve a system of 
quadratic equations (anyway, the computation time is a few seconds). The results are shown 
in Tables [l] [2] [3j [4] (until reaching a relative error of order 10~^ or less) and Figure [4j 



n 


•5 71 


cxpp s„ 


Rel. 


error 





(0,0,0,1) 


(0, 10.65556, 7r/2, 0.821858) 


2.1 




1 


(0,-5.928809,0,0.835241) 


(0, 5.981945, 7r/2, 1.421563) 


5.9 • 


10-1 


2 


(0,-7.220423,0,0.604050) 


(0, 4.050859, 7r/2, 1.568339) 


1.4 • 


10-2 


3 


(0,-7.247942,0,0.597185) 


(0, 4.000075, 7r/2, 1.570793) 


2.1 ■ 


10-5 


4 


(0,-7.247982,0,0.597175) 


(0, 3.9999998, 7r/2, 1.570796) 


6.2 ■ 


10-8 



Table 1 
expp 

the relative errors of each coordinate between exp^ s. 
left). 



: Example 4.1 Successive iterations of the algorithm for computing s. The event 
approximates the intersection event ^s- The relative error corresponds to the sum of 

and gs = (0, 4, 7r/2, 7r/2) (see Figure [i] 
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n 




expp s„ 


Rel. 


error 





(0,0,0,1) 


(0, 10.65556, 7r/2, 0.821858) 


2.1 




1 


(0,-6.940805,0,0.666507) 


(0,4.531059,71/2, 1.538718) 


1.5 • 


10-1 


2 


(0,-7.247475,0,0.597303) 


(0, 4.000944, 7r/2, 1.570752) 


2.6 ■ 


10-4 


3 


(0,-7.247982,0,0.597175) 


(0, 3.99999991, ir/2, 1.570796) 


2.9 ■ 


10-8 



Table 2: Example |4.1[ Analogous to Table [T] but using the second order approximation of 
Remark 13.31 



n 


Wr 


I 


Cxpp Wn 


Rel. 


error 





(- 


9.237604,0,0, 1) 


(- 


8.991730, 10.245645, n/2, 0.844573) 


2.1 




1 


(- 


6.445192, -3.916650, 0, 0.408921) 


(- 


7.241158, 4.720295, tt/2, 0.723016) 


3.9 ■ 


10-1 


2 


(- 


6.473183, -4.365314, 0, 0.306661) 


(- 


7.571348, 4.020019, 7r/2, 0.627855) 


1.4 ■ 


10-2 


3 


(- 


6.472659, -4.3769972, 0, 0.303040) 


(- 


7.580813, 4.000015, 7r/2, 0.623213) 


1.4 ■ 


10-^ 


4 


(- 


6.472658, -4.377011, 0, 0.303036) 


(- 


7.580825, 3.999991, 7r/2, 0.623207) 


2.2 ■ 


10-6 



Table 3: Example |4.1| Successive iterations of the algorithm for computing w. The event 
exppM;„ approximates the intersection event qi. The relative error corresponds to the sum 
of the relative errors of each coordinate between exp^Wn and qi (see Figure |4] right). In this 
case, since the exact t and ip coordinates of qe are unknown, they have been assumed to be 
the t and ip coordinates of exp 104 ; the r and 6 coordinates of qi are 4 and tt/2 respectively. 



Example 4.2 The Kerr metric in Boyer-Lindquist coordinates {t, r, 9, (p} is given by the line 
element 

ds^ = -df + ^dr^ + p^de^ + (r^ + a^) sin^ Od^^ +'^r{dt-a sin^ 9d^) \ (17) 



where 



A " ■ V ■ / p 



A - 2'mr + a^, + cos^ ( 



This metric describes the exterior gravitational field of a rotating mass m with specific angular 
momentum a = J /m^ where J is the total angular momentum of the gravitational source 
(see 19 for restrictions on the coordinates). From now on we are going to suppose that 



m = 1 and a = 1/2. 

Analogously to Example 4.1 given a stationary observer at ro = 8, 6*0 = 7r/2, (p^ = 0, 
and a test particle with equatorial circular geodesic orbit with ri — A, 9i — 7r/2, Lpi — 7r/2 at 
t = Q, we are going to check the algorithm for computing s and w aX p = (0, 8, 7r/2, 0) using 
an initial spacelike direction sp = (~l/3, 0, 0, 1) (orthogonal to the 4- velocity of the observer 
at p), and using an initial past-pointing lightlikc direction wo = (—1/3 — ^3091/6, 0, 0, 1), 
respectively. The method works as well as in Example |4.1[ moreover, in the lightlike case, 
the linear estimation gives similar results to those applying the second order approximation 



proposed in Remark 3.3 The results are shown in Tables [sj [6] [t] and [s] (until reaching a 



relative error of order 10 or less). 



n 


It!, 


1 


expp Wn 


Rel. 


error 





(- 


-9.237604,0,0,1) 


(-8.991730, 10.245645, 7r/2, 0.844573) 


2.1 




1 


(- 


6.466102, -4.131156, 0, 0.366627) 


(-7.389968, 4.398864, 7r/2, 0.691871) 


2.3 


10-1 


2 


(- 


6.472884, -4.372292, 0, 0.304508) 


(-7.576998, 4.008086, 7r/2, 0.625106) 


5.6 


10-3 


3 


(- 


6.472658, -4.377008, 0, 0.303037) 


(-7.580822, 3.999996, 7r/2, 0.623209) 


1.0 


10-6 



Table 4: Example |4.1| Analogous to Table [3) but using a second order approximation 
analogous to that of Remark |3.3[ 
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n 






expp Sn 


Rcl. 


error 





(- 


0.333333,0,0,1) 


(- 


0.244901, 10.664776, 7r/2, 0.822363) 


2.9 




1 


(- 


0.266769, -5.872871, 0, 0.800307) 


(- 


0.706191, 5.752058, 7r/2, 1.386445) 


8.4 


10-1 


2 


(- 


0.197485, -6.878255, 0, 0.592454) 


(- 


1.084179, 4.026634, 7r/2, 1.450633) 


1.5 


10-2 


3 


(- 


0.196280, -6.889558, 0, 0.588841) 


(- 


-1.092323, 4.000189, 7r/2, 1.450052) 


1.0 


10-4 


4 


(- 


0.196272, -6.889639, 0, 0.588815) 


(- 


-1.092382, 3.9999995, 7r/2, 1.450047) 


1.2 


10-^ 



Table 5: Example 4.2 Successive iterations of the algorithm for computing s. The event 
expp s„ approximates the intersection event qs- The relative error corresponds to the sum of 
the relative errors of each coordinate between exp^ s„ and — (0, 4, 7r/2, 7r/2) (see Figure 
|4]left). In this case, since the exact t and coordinates of q^ are unknown, they have been 
assumed to be the t and coordinates of exp^ 54; the r and 6 coordinates of qs are 4 and 
7r/2 respectively. 



71 






Cxpp Sn 


Rcl. 


error 





{- 


0.333333,0,0,1) 


(-0.244901, 10.664776, 7r/2, 0.822363) 


2.9 




1 


(- 


0.215246, -6.683244, 0, 0.645738) 


(-0.968833, 4.433499, 7r/2, 1.449471) 


2.2 • 


10-1 


2 


(- 


0.196243, -6.889900, 0, 0.588730) 


(-1.092573, 3.999382, 7r/2, 1.450033) 


3.4 • 


10-4 


3 


(- 


0.196272, -6.889638, 0, 0.588815) 


(-1.092382, 3.99999995, 7r/2, 1.450047) 


1.1 ■ 


10-8 



Table 6: Example |4.2[ Analogous to Table [5j but using the second order approximation of 
Remark 13.31 



Next, we are going to give some examples in Schwarzschild and Kerr spacetimes about 
computing the relative velocities along an observer. In the figures, we make "retarded com- 
parisons" (see [To]) of their moduli, i.e. we plot the modulus of a relative velocity as function 
of ts (kinematic and Fermi) or ti (spectroscopic and astrometric) . All the numerical data has 
been computed with a relative error less than 10~^. 



Example 4.3 In the Schwarzschild metric (16), let us consider the stationary observer and 



the test particle of Example |4.1[ but now we are going to suppose that it has tpi — at t = 
(i.e. the observer, the test particle and the singularity r = are aligned at t = 0). This 
problem has not been previously studied analytically due to its complexity. 

Note that in this case, we do not work in a convex normal neighborhood and so there is not 
a unique kinematic, Fermi, spectroscopic and astrometric relative velocities, depending on 
different geodesies joining the observer and the test particle (i.e. different choices of ip in the 
spacelike case, or A in the lightlike case, see Figure [T]). For example, if the test particle is at 
spatial coordinates ri = 4, 0i = 7r/2, (pi = 0, the observer and it can be joined by geodesies, 
■0 or A, giving whole turns around the black hole or not. In Figure [5] there are represented the 
moduli of the corresponding relative velocities in the case of equatorial geodesies joining the 
observer and the test particle following the convention that ipi also indicates the number of 
turns (and their direction) around the black hole: for example, for pi ~ 0, the geodesic goes 
directly from the observer to the test particle; on the other hand, for ipi = 27r the geodesic 
gives an equatorial whole turn counter-clockwise around the black hole before arriving at the 



n 


w. 


1 


expp w„ 


Rel. 


error 





(- 


9.599460,0,0,1) 


(- 


9.288415, 10.147354, tt/2, 0.854357) 


2.0 




1 


(- 


6.805715, -4.000178, 0, 0.432700) 


(- 


7.945461, 4.553253, 7r/2, 0.747437) 


2.9 ■ 


10-1 


2 


(- 


6.764920, -4.337490, 0, 0.356268) 


(- 


-8.176669, 4.005345, tt/2, 0.667962) 


3.4- 


10-3 


3 


(- 


6.764068, -4.340567, 0, 0.355373) 


(- 


-8.178553, 3.999982, 7r/2, 0.666767) 


4.4 ■ 


lo-'' 



Table 7: Example 4.2 Successive iterations of the algorithm for computing w. The event 
exppM;„ approximates the intersection event qi. The relative error corresponds to the sum 
of the relative errors of each coordinate between exp^Wn and qg (see Figure |4] right). In this 
case, since the exact t and ip coordinates of q^ are unknown, they have been assumed to be 
the t and ip coordinates of exp^ w^] the r and 6 coordinates of qi are 4 and tt/2 respectively. 
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n 




expp Wn 


Rcl. 


error 





(-9.599460,0,0,1) 


(-9.288415, 10.147354, 7r/2, 0.854357) 


2.0 




1 


(-6.747007, -4.393588, 0, 0.339049) 


(-8.209192, 3.906149, 7r/2, 0.643900) 


6.1 


10-2 


2 


(-6.764157, -4.340248, 0, 0.355466) 


(-8.178359, 4.000538, 7r/2, 0.666892) 


3.4 


10-4 


3 


(-6.764068, -4.340567, 0, 0.355373) 


(-8.178554, 3.999982, tt/2, 0.666767) 


4.4 


10'^ 



Table 8: Example |4.2| Analogous to Table [7j but using a second order approximation 
analogous to that of Remark |3.3[ In this case, the result for n = 3 is very similar to the case 
of the standard algorithm (see Table [7| . 




Figure 5: Retarded comparison of the moduli of the kinematic, Fermi, spectroscopic and 
astrometric relative velocities of a test particle with equatorial circular geodesic orbit with 
radius ri = 4, 6*1 = 7r/2 and tpi = &t t = 0, with respect to a stationary observer at ro = 8, 
^0 = 7''/2 and ipa — 0, in the Schwarzschild metric with m = 1. The vertical lines correspond 
to different values of ipi (the (/7-coordinate of the test particle at or qi), and the geodesies 
joining the observer and the test particle have been restricted to be equatorial. 



test particle. 

At t = 0, we have taken initial vectors sq — (0,-1,0,0) and wq ~ (—4/3,-1,0,0) (for 
being a past-pointing lightlike vector at p = (0, 8, 7r/2, 0)), but this is not important because 
in this case the algorithm converges quickly also for non-nearby initial vectors and so it is not 
necessary to apply the second order approximation of Remark |3.3[ The computations have 
been done using a time step (in the observer) of At — 0.25, and the number of iterations 
needed at each time (for reaching the desired relative error 10"®) is at most n — i. So, the 
linear algorithm works very well in this case. 

Considering stationary observers with different radial coordinate = 4 and tq = 3 (see 
Figures [g] and [t] respectively), we observe that ||wkin|| remains constant and equal to ^1/2. 
This numerical result has motivated a work [20] where this property is theoretically proved 
in general: in the Schwarzschild metric, the modulus of the kinematic relative velocity of a 
test particle with circular geodesic orbit at radius ri > 3m with respect to any stationary 
observer is constant and equal to . / — . 

Moreover, it can be checked that ||uspoc|| tends to 1 when ti ±oo, and using expression 
([3]), we can compute the corresponding frequency shift, as it is seen in Figure p| 
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Figure 6: Analogous to Figure [5] but taking rg — 4. 



1^1 = — 27r t/'i — — TT i^i = i/^i — 7r = 27r 
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Figure 7: Analogous to Figure [Sj but taking rg = 3. 
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Figure 8: Frequency shifts corresponding with the spectroscopic relative velocities computed 
in Example 4.3 (see Figures [5)[6]andl7|. 



Example 4.4 In the Kerr metric ( |17[ ), let us consider the stationary observer and the test 
particle of Example |4.2[ but supposing that it has ip\ = at t = 0, as in Example |4.3[ Anal- 
ogously to this example, in Figure [9] there are represented the moduli of the corresponding 
relative velocities. Comparing with the analogous problem in the Schwarzschild spacetime 
studied in Example 4.3 and Figurcjsj it can be observed that |lwkin|| does not remain constant. 
Moreover, it also draws attention to the fact that ||fspoc|| does not tend to 1 when tg — > — oo; 
in fact, it can be checked that it is decreasing and tends to a value « 0.072. 

Finally, in Figure [T0| it is shown the frequency shift of the test particle with respect to the 
observer, compared with the frequency shift of the analogous problem in the Schwarzschild 
spacetime. Of note is the fact that, in the Kerr spacetime, the shift is greater than 1 for a 
sufficiently negative Lp\\ in fact, it tends to 1.053 approximately when (^i — > —oo. Hence, in 
this case, an approaching?] test particle has redshift instead of blueshift. 



5 Final remarks 

First, we have generalized the concepts of the relative velocities introduced in [7 to non- 



convex normal neighborhoods (see Section 2.1), focusing on what is minimally necessary to 



make sense of the corresponding definitions. As a result, we can now apply this theory of 
relative velocities to a wide range of scenarios, including those with gravitational lensing or 
caustics. 

Then, we have developed an algorithm for computing relative velocities of a test particle 
with respect to an observer, based on computing the Fermi and observational coordinates of 
the test particle; hence, this method can be applied only for this purpose. With respect to 



Fermi coordinates, the objective of the paper 21 is also to find the Fermi coordinates of an 



event in a general spacetime, calculating the general transformation formulas from arbitrary 
coordinates to Fermi coordinates, and vice versa. But the methods are not the same: 

• In [ 2Tj , the Fermi coordinates are given in the form of Taylor expansions. Hence, for 
increasing accuracy, we have to add terms of higher order, and it could be very complex 
to achieve a high accuracy for distant events. 



■^Taking into account the affine distance, also known as lightlike distance, defined as ||sobsll (see 22 ). 
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Figure 9: Retarded comparison of the moduli of the kinematic, Fermi, spectroscopic and 
astrometric relative velocities of a test particle with equatorial circular geodesic orbit with 
radius ri — A, 9i = 7r/2 and ipi — at t — 0, with respect to a stationary observer at vq = 8, 
Oa = 7r/2, (fio = 0, in the Kerr metric with m = 1, a = 0.5. The vertical lines correspond 
to different values of ipi (the (/3-coordinate of the test particle at Qs or qi), and the geodesies 
joining the observer and the test particle have been restricted to be equatorial. 




Figure 10: Frequency shifts corresponding with the spectroscopic relative velocities computed 
in Example 4.4 (see Figure [9]) and in the analogous problem in the Schwarzschild spacetime 



(see Figure [5| and Figure [s] with rp = 8). The frequency shifts are plotted as functions of (pi 
(the (/5-coordinate of the test particle at q^) in order to make a fair comparison. 
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In this work, we use an iterative algorithm. So, for increasing accuracy, we have to 



make more iterations (provided there is convergence, see Remark 3.4), and this does 
not imply additional difficulty. 

Moreover, the algorithm proposed in this work is also valid in non-convex normal neighbor- 



hoods where there is not uniqueness of Fermi coordinates (see Section 3.4): depending on 
the initial vector sq you can get different relative positions s of the same test particle at 
a given observer's time. Nevertheless, the method given in [21| is a powerful tool in the 
theoretical study of Fermi coordinates, while the algorithm introduced here is designed to 
be implemented on computers. So, they are complementary and, for example, we can use 



the Taylor expansions given in 21 for choosing the initial vector sq in the first execution. 



Moreover, if the algorithm does not converge (see Remark 3.4), then the Taylor expansions 
are a valuable alternative. 

Finally, it is important to remark that there exist some fundamental applications derived 
from the computation of the relative velocities, e.g. we can find the frequency shift and 
the light aberration effect (see |22:) from the spectroscopic relative velocity, or measure the 
expansion of space from the kinematic and Fermi relative velocities of comoving observers of 



FLRW spacetimes (see 11 12^). 



A Alternative computation of Fermi and astrometric 
relative velocities 



In 23 Proposition 3.3] it is proved that, working in a convex normal neighborhood, vpeimi 
and Vast can be computed in terms of p, Qs, qe, u, (VuC/)^, u'^, u'^, s, Sobs, and hence we do 
not need to know 5* or S'obs around p, contrary to what is expected from Q and This 
computation is done using a coordinate system ( ) , obtaining 

tTermi = + T a2 + a-^ - g [s, {VuU)^^ u, (18) 
where the coordinates of vectors ai, a2, a^, are given by 

a't-.^dJtiPK ; a^=a,logj^(g,K'' ; a[; :^ T^pK , 

with 

/s :=log(_,gs), (19) 
where log(p, g) denotes log^q (or equivalently exp~^ q), and 



g (^s, (V(7f/)p) +g{ai+ 03, u) 



g(a2,u) 

On the other hand 

Wast = 04 + 'r'a5 + (3(04 + 'r'a5,") +5 (logp(jf,-u)) u 

+9 (logp qi, u)u + aG- g (^Sobs, (Vc/t/)^) u, (20) 

where the coordinates of vectors 04, 05, ag are given by 

4:=d^fjf{pK ; a^=a.log|^(g,)«r ; a^, -.^ r^PK s^^^, 

with 

/£:= log (_,?£), (21) 



18 



and 



5 (sobs, (Vc/C/)p) + g (ai + a6,u) + g (\ogpqt, g iu,u) ' 



9 (04, w) 

Note that ii is given in terms of p, u and {WuU)^: 



Expressions (18) and (20 1 of WFcrmi and Wast are not explicitly shown in [23', but they can be 



deduced from the proof of Proposition 3.3. 

Nevertheless, numerically it is difficult to compute the vectors oi and 04 with high accu- 



racy, concretely the derivatives of /s (19) and Ji (21). For example, if we want to compute 



f^/i/si we are supposed to know qs and the vector logp^g, i.e. the relative position s, that can 
be estimated by means of the algorithm exposed in Section |3.2[ Then, given a small e > 0, 
we launch a geodesic from the event p' with coordinates p^ + eJ^ (i.e. the same coordinates 
as p but the /U-th coordinate is p^^ + e) and initial tangent vector s (actually, the vector in 
TpiM. with the same coordinates as s). Since e is small, this geodesic is "close"[^to qs, and 
so there is an event ggeo of the geodesic "close" to qs- So, assuming an afhne structure, we 
have to parallel transport the vector qs — (?geo from q^eo to p' along the geodesic, and add 
this vector to the current initial vector s, obtaining a new initial vector whose corresponding 
geodesic will be "closer" to qs- Repeating this process, we can estimate the initial vector of 
the geodesic passing through qs with the desired accuracy and then, we can evaluate d^fs 
comparing this initial vector with s. 

Analogously, for computing the derivatives of fi, we are supposed to know q^ and the 
vector log^g^, whose projection onto is the observed relative position Sobs, that can be 
estimated by means of the algorithm exposed in Section [3.3[ 

Concluding, this method let us find i^Formi and v^st computing S and S'obs only at p, but 
the original method based on definitions (|4| and ([s]) (in which 5' and S'obs are computed 
around p) is obviously faster and more accurate because it requires a far fewer number of 
operations. Moreover, if we do not work in a convex normal neighborhood, expressions ( 18 ) 



and (20) are not strictly valid because the vectors ai, 02, 04, and 05 are not well-defined in 
general. 
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